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ABSTRACT 

Three  probability  models  for  defense  against  mass  missile  attacks  are 
developed.  Each  model  corresponds  to  a  different  level  of  information.  Different 
strategies  are  made  possible  by  the  information  available.  The  first  strategy,  Invisible 
Kill,  assigns  equal  times  to  the  destruction  of  each  missile  that  is  selected  for 
attention.  The  second,  a  Threshold  strategy,  assigns  a  maximum  threshold  time  to 
the  same  task.  The  last,  based  on  Perfect  Task  Information,  engages  missiles  in 
ascending  order  of  their  destruction  times;  the  order  being  assumed  known  in 
advance.  The  Expectation  of  the  number  of  missiles  killed  is  the  Measure  of 
Effectiveness  used  to  evaluate  the  models.  Numerical  results  are  analyzed  through 
simulation.  Different  strategies  are  compared  showing  the  effects  of  information  in 
defense  effectiveness. 
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THESIS  DISCLAIMER 

The  reader  is  cautioned  that  computer  programs  developed  in  this  research  may 
not  have  been  exercised  for  all  cases  of  interest.  While  every  effort  has  been  made, 
within  the  time  available,  to  ensure  that  the  programs  are  free  of  computational  and 
logic  errors,  they  cannot  be  considered  validated.  Any  application  of  these  programs 
without  additional  verification  is  at  the  risk  of  the  user. 
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I.    INTRODUCTION 

Information  is  known  to  be  an  important  factor  in  the  determination  of  the 
outcome  of  military  combat.  Consequently,  the  addition  of  information-gathering  and 
interpretation  assets  must  be  viewed  as  in  legitimate  competition  with  the  acquisition 
of  weapons.  In  fact,  modern  "smart"  weapons  possess  within  themselves  a 
combination  of  brawn:  destructive  power  (e.g.  high  explosive  warhead);  and  brain: 
navigation,  identification,  guidance  and  homing  power  (e.g.  on-board  sensors  and 
logic).  External  C2  systems  are  needed  to  guide  overall  decision-making:  on  the 
defender's  side  to  detect  attacks  and  to  schedule  defense  response;  on  the  attacker's 
side  to  first  soften  or  confuse  (EW  measures,  deception)  the  defender  prior  to  a 
focused  destructive  action  against  sensitive  assets. 

This  thesis  seeks  to  investigate  the  effects  and  value  of  information  availability 
in  a  defense  against  a  single  mass  (simultaneous)  missile  attack.  The  analysis  is 
carried  out  from  the  point  of  view  of  a  defender  operating  within  constraints  of  a 
fixed  time  window,  presumably  until  destruction  of  the  attacking  missiles  is  no 
longer  possible  at  the  level  considered;  "leaking"  missiles  may  be  dealt  with  at  an 
ensuing  defense  level.  For  simplicity  we  postulate  that  the  objective  of  the  defender 
is  to  maximize  the  expected  number  of  missiles  destroyed  at  the  level  under 
consideration.  The  defender  is  using  some  weapon  system  (interceptor  missiles, 
irradiation  weapons,...).  This  generic  weapon  system  is  characterized  by  its 
probability  of  acquiring  a  target,  acquisition  or  setup  time,  and  rate  of  kill.  These 
parameters  are  supposedly  design  features  intrinsic  to  the  system.  It  is  possible  that 
trade-offs  between  the  weapon  system  characteristics  (parameters)  exist,  in  which 


case  our  models  may  provide  some  insight  in  deciding  the  desirable  parameter 
configuration. 

The  three  models  we  develop  are  of  a  probabilistic  nature,  each  of  them 
corresponding  to  a  different  information  level  available  to  the  defender.  A  higher 
information  level  should  enable  the  user  to  enhance  his  strategy  and  consequently 
achieve  a  gain  in  defense  effectiveness.  The  general  methodology  followed  is  to  start 
with  the  basic  assumptions  of  the  model  and  strategy,  and  work  our  way  analytically 
to  explicit  expressions  for  the  chosen  measures  of  effectiveness.  Numerical  results 
of  the  probability  models  are  checked  against  those  of  discrete-event  simulations. 

It  must  be  noted  that,  in  this  context,  the  term  "information"  is  used  in  a 
specific  but  informal  sense  appropriate  to  the  present  application.  There  is  no  direct 
connection  with  such  classical  concepts  as  Shannon  information  or  Fisher 
information.  Finally,  note  that  the  present  problem  may  appear  in  other  guises,  in 
which  a  number  of  tasks  must  be  performed  under  time  constraints;  e.g.,  in  repair 
and  maintenance  or  flexible  manufacturing. 


II.    MASS  ATTACK  WITH  INVISIBLE  KILL 

A.  INTRODUCTION 

This  chapter  deals  with  the  following  situation:  A  defender  is  under  a  mass  or 
simultaneous  attack  from  a  possibly  large  number  of  incoming  missiles  (a  "cluster"). 
Available  to  the  defender  is  an  unspecified  outer-layer  weapon  system,  Dl,  which 
enables  him  to  shoot  at  least  some  of  the  missiles.  Each  missile,  upon  being  detected 
must  be  acquired,  and  only  then  can  one  attempt  to  shoot  it  down.  The  weapon 
system  only  allows  the  defender  to  engage  one  target  at  a  time,  so  he  faces  the 
problem  of  how  to  optimally  assign  time  to  each  missile  so  as  to  maximize  some 
measure  of  the  number  of  missiles  destroyed.  It  is  understood  that  one  is  working 
in  a  restricted  time  window:  if  not  destroyed,  the  missiles  will  proceed  towards  their 
targets,  passing  through  a  secondary  defence  level.  The  missiles  that  are  not  killed 
by  Dl  are  said  to  leak  to  the  next  layer.  The  general  objective  is  to  minimize 
leakage. 

B.  ASSUMPTIONS 

In  the  probability  model  the  following  simplifying  assumptions  will  be  used: 

•  The  number  of  attacking  missiles  is  known. 

•  The  probability  of  acquiring  each  target  is  known  and  fixed;  range  dependence 
is  ignored. 

•  Given  that  the  missile  is  acquired,  the  probability  distribution  of  the  time  to 
destroy  it  is  known. 

•  It  is  not  known  to  the  defender  whether  a  given  missile  was  successfully 
acquired  or,  if  acquired  was  actually  killed.  This  is  the  minimal  information 
situation;  we  call  it  Invisible  Kill. 


•    There  is  a  fixed  setup  time  expended  each  time  a  new  attempt  is  made  to 
acquire  another  missile. 

Finally,  the  strategy  of  the  defender  will  be  to  allocate  time  evenly  between  those 
missiles  selected  for  attention  by  D 1 .  Not  all  attacking  missiles  will  be  chosen  under 
optimal  conditions;  some  will  be  allowed  to  leak,  presumably  with  warning 
transmitted  to  the  second  (inner)  defence  layer. 

C.      MODEL 

The  following  notation  will  be  used:  Let  n  be  the  number  of  incoming  missiles. 
The  probability  of  acquiring  a  target  is  p  ,  and  q  -  \-p  the  probability  of 
acquisition  failure.  Denote  by  X  the  random  time  to  destroy  a  missile  which  has 
been  acquired.  T  is  the  time  available  for  destroying  the  targets  -  the  duration  of 
the  window  of  opportunity.  There  is  a  fixed  time  necessary  for  each  acquisition 
which  will  be  denoted  by  w. 

It  is  assumed  that  the  distribution  of  X  is  known.  Call  its  CDF  Fx(x) ,  where: 

X  ~  Fx(x)  =  Prob  {  X  <>  x  \  Target  was  acquired  )  .  (2.1) 

Suppose  one  chooses  to  engage  a  number,  m,  of  the  incoming 
missiles;  m  <>  n  .  Recalling  the  declared  strategy  of  assigning  equal  time  to  each 
one,  then  T  -  mw  +  mx;  here  x  is  the  time  allowed  to  destroy  a  single  missile. 
It  follows  that: 


x  =  -  -  w  .  (2.2) 

Since  each  target  must  be  destroyed  in  less  than  time  x,  Prob  { destroying  a 
target}  =  Prob  (a  target  is  acquired  upon  attempt}  *  Prob  (a  target  is  destroyed  in  the 
time  available:    max(0,  T/m-w)     };  i.e.: 

T 
Prob  {destroying  a  target)  =  p  .  Fx(  —   -  w  )  .  (2.3) 

m 

Note  that  since  X  2  0  with  probability  one,  automatically    Fx(x)  -  0   if  x  s  0  . 

If  it  is  assumed  that  each  attempt  to  acquire  and  destroy  a  target  is  an 
independent  Bernoulli  trial,  then  the  number  of  targets  destroyed  out 
of  m,  N(T),  will  be  distributed 

N(T)  ~  Binomialim  ,pFx(—-w))y  (2.4) 

m 

and  the  expected  number  of  targets  destroyed  will  be 

E[N(T)]  =  mpFx(—-w).  (2.5) 

m 

The  current  problem  is  to  find  the  optimal  value  of  m ,  i.e.,  one  that  maximizes  this 
expected  value.  Other  criteria  and  controls  are  of  course  possible.  Note  that  if  the 
targets  are  known  to  be  different,  then  a  strategy  that  schedules  different  times  to 
each  type  would  be  reasonable.  An  extension  of  this  formulation  might  include  a 
preliminary  classification  of  target  type  and  subsequent  allocation  of  effort  to  the 
most  profitable  types  first. 


D.      MODEL  DEVELOPMENT 

The  optimal  value  of  m  will  now  be  derived:  namely  the  number  of  missiles 
the  defender  should  engage  in  order  to  destroy  the  maximal  expected  number  of 
them  at  his  particular  defense  layer. 

The  expected  number  of  destroyed  or  killed  missiles  is  given  by  (2.5).  We 

wish  to  find  the  maximum  of  this  function.  Set  the  derivative   *■   -     j    =  0 , 

dm 

implying  that  if  the  optimal  value,  m  =  m  °,  is  interior,  0  <  m  °  <  n,  then 

Fx(z-w)  =  Zfx(z-w),  (2.6) 

where  the  parameter  z  =  T/ m°.   In  principle  there  can  be  several  such  stationary 
points,  and  furthermore  the  optimal  number  can  be  a  boundary  value:   m°  -n. 

In  particular  suppose  X  -  Exponential  (A) ,  where  the  rate  parameter  may  be 
interpreted  as  the  number  of  acquired  missiles  the  weapon  kills  per  unit  time.  Then: 

Fx  (x)  =  1  -  e  ~Xx,     x  >  0  ;  (2-7) 

The  expected  number  of  destroyed  missiles  becomes: 

-*(z-w)  (2.8) 

E[N(T)]  =  mp[l   -  e      m      ].  v      ' 

To  maximize  this  function  of  m,  take  the  derivative  and  set  it  equal  to  zero. 
After  some  algebra,  get: 

ln(m)  -  ln(m  +  AD  +  —  -  Aw  =  0  ;  (2.9) 

m 

An  approximate  explicit  solution  for  m  * ,  a  stationary  point,  even  if  somewhat 
inaccurate,  would  be  desirable  to  show  its  dependency  on  the  various  parameters. 
Start  by  rewriting  (2.9),  and  using  a  Taylor  expansion: 


XT 

—  -  kw  =    ln(m+A.r)  -  lnm 

m 

m  m         2m2 

solving  for  m,  yields  the  approximate  explicit  result: 


(2.10) 


m 


A_  (2.11) 


This  number  will  be  a  positive  real.  We  want  an  integer  solution;  therefore  the 
greatest  integer  not  exceeding    m  * ,  [m*],    and  the  next  lowest  integer,  <m  '  >, 
must  be  checked  to  see  which  one  has  the  largest  objective  function  (2.8).  Call  this 
number  m**,  then  the  approximate  optimal    m°  will  be: 

m°  =  min(rt,m**).  <2-12) 

Accordingly  the  corresponding  optimal  value  of  the  expected  number  of 
missiles  killed  is,  again,  approximately: 

-H-0-w)  (2.13) 

E°[N(T)]  =  m°p  [l-e    u       ']. 

Note  that  the  second  derivative  of  E[N(T)]  is 

d2E[N(T))   _        (AT?    -*(j-»).  (214) 

dm2  m3 

therefore 


^E[N(T)]<QVm  (215) 

dm2 

The  function  E[N(T)]  is  concave  throughout  the  domain  of  m.  In  this 
case  m  °,  if  obtained  as  described  above,  will  be  the  approximate  unique  global 
maximizer. 

E.      SIMULATION 

To  verify  the  analytical  results  obtained  thus  far,  a  terminating  simulation 
[Ref.  1]  was  performed.  The  measure  of  effectiveness  estimated  is  the  Expected 
Number  of  Missiles  Killed  per  Attack.  The  attacks  were  replicated  independently, 
15000  times,  with  several  combinations  of  input  parameters.  The  FORTRAN  codes 
pertaining  to  the  Event  routines  are  displayed  in  Appendix  A.  For  Timing  and 
Calendar  manipulation  routines  the  SIMUTIL  package  [Ref.  2]  was  employed,  as 
was  the  LLRANDOMII  Random  Number  Generator  [Ref.  3]. 

The  results  of  the  simulation  were  written  to  CMS  files,  and  used  to  generate 
several  plots  of  response  surfaces.  One  of  these  plots  is  shown  in  Figure  2.1,  with 
E[Missiles  Killed]  as  the  response  variable. 

For  example,  if  one  takes  the  time  window  to  be  700,  20  attacking  missiles, 
the  rate  of  kill  0.01,  the  acquisition  time  10  and  the  probability  of  acquisition  0.85; 
then  solving  (2.11)  yields  m  *  =  15.652476.  Checking  m*  =  15,  16  with  (2.13) 
reveals  that  m  *  *  -  15 ;  therefore  m  °  =  min(  15, 20)  =  15 .  It  can  be  observed  from 
Figure  2.1,  with  these  values,  that  the  approximate  analytical  optimal  seems  to  be 
in  good  accordance  with  the  simulation  results.  Note  that  obtaining  the  precisely 
correct  value  of  m  is  not  important;  the  value  of  the  expected  number  killed  is  not 
very  sensitive  to  that  value. 
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Figure  2.1   Simulation  Results 


F.      USE  OF  DECOYS  BY  ATTACKER 

It  might  be  of  some  practical  interest  to  discuss  the  influence  of  the  use  of 
decoys,  in  the  results  of  the  previous  strategy  for  Invisible  Kill.  By  assumption  the 
decoys  will  be  randomly  mixed  with  the  missiles  within  the  cluster.  As  information 
available  is  minimal,  the  decoys  will  be  undistinguishable  from  the  missiles  for  the 
defender;  as  a  result  the  defender  employs  the  same  strategy  regardless  of  their 
presence. 

Let  r  be  the  fraction  of  targets  which  are  true  missiles,  i.e., 


number  of  missiles  (2  16) 

number  of  missiles  +  number  of  decoys 

Using  the  current  independence  assumptions,  the  number  of  missiles  destroyed 
in  m   trials  will  be  distributed  Binomial  and  its  mean  will  be: 

T 

E[missiles  destroyed]  =  mprFx{—   -  w) .  (2.17) 

m 

If  the  destruction  time  is  again  taken  to  be  Exponential,  the  same  steps  already 
taken  show  that  the  optimal  number  of  targets  to  engage  does  not  vary  with  r.  The 
expected  number  of  missiles  destroyed  using  the  optimal  strategy  decreases,  being 
multiplied  by  a  factor  of  r.  This  result  is  rather  intuitive.  Since  there  is  no  way  the 
defender  can  assess  the  composition  of  the  attacking  target  cluster,  defense 
effectiveness  depends  directly  on  the  fraction  of  decoys.  Eventually  an  attack  with 
large  numbers  of  decoys  will  succeed  in  diluting  the  defender's  effort,  thus 
achieving  a  high  rate  of  penetration  to  a  subsequent  defense  layer. 

It  is  interesting  to  conjecture  the  payoff  from  a  subsystem  capable  of 
classifying  attacking  entities  as  either  missiles  or  decoys  with  some  probability 
("skiH").  It  may  be  that  tradeoff  between  some  of  the  other  parameters, 
e.g.,  T,  w,  or  X,  and  classification  skill  can  be  beneficial  for  the  overall  system 
effectiveness.  This  tradeoff  is  susceptible  to  quantitative  study  conducted  in  the 
mode  described  throughout  this  thesis. 
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HI.    VISIBLE  KILL  -  A  THRESHOLD  STRATEGY 

A.  INTRODUCTION 

The  former  chapter  dealt  with  a  situation  in  which  there  is  no  direct 
information  pertaining  to  when,  or  if,  a  given  target  is  destroyed.  The  previous 
assumptions  are  next  modified  in  a  way  such  that  the  defender  recognizes  the  instant 
when  a  target  is  destroyed  (VISIBLE  KILL),  if  such  an  event  occurs.  Upon 
destruction  of  a  target  the  defender  immediately  attempts  to  acquire  another  (if  any 
are  left).  It  should  be  noted  that  acquisition  is  still  'not  visible',  i.e.  after  the  fixed 
acquisition  time  w,  the  defender  does  not  know  if  the  acquisition  was  a  success. 
Thus  Dl  may  have  to  wait  indefinitely  for  news  that  a  kill  has  occurred,  which  leads 
to  consideration  of  a  time  threshold  strategy  or  policy. 

B.  THRESHOLD  MODEL 

A  threshold  policy  is  adopted:  the  defender  will  allow  a  maximum  threshold 
time  ,  x,  to  be  allocated  to  destruction  of  a  target.  If,  after  time  x  elapses,  the 
target  is  not  observed  to  be  destroyed  the  defender  will  go  back  to  acquisition  mode. 
Such  a  policy  is  easy  to  implement  and  can  be  justified  by  dynamic  programming 
under  certain  circumstances. 

Define  Z(x),  the  total  random  time  to  kill  one  target  with  a  threshold  equal 
to   i.   Then,  from  the  assumptions  made,  the  following  recursion  must  hold: 
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Z(t)|X=x  =  < 


w  +  x  if  x  <>x 

(3.1) 

> 

w  +  x  +  Z'(x)     if  x>  x 


where  X  is  the  random  time  to  destroy  an  acquired  target,  in  a  continuous 
engagement.  Z  (x)  is  a  random  variable  having  the  same  distribution  as  Z(t);  it 
is  the  time  to  shoot  down  a  target  if  the  first  attempt  fails  and  the  process  of 
acquisition  must  start  over.  The  CDF  of  X  is  known  to  be  Fx(x). 

The  upper  branch  of  (3.1)  occurs  when  acquisition  is  successful  and  the  target 
is  shot  within  the  prescribed  threshold.  This  has  probability  pdFx(x).  The  lower 
branch  occurs  if  either  acquisition  fails,  or  acquisition  succeeds  but  the  missile  is  not 
shot  in  the  allowed  threshold  time.  This  has  probability  q  +  pFx(x).  Thus  we 
obtain  the  unconditional  expectation: 


£[Z(x)]  =  °- 


w  + fpxdFx(x)  +x[q +pFx(x)]  (3  2) 


pFJx) 


where  Fx(x)  =  1   -  Fx(x) .  Integration  by  parts  in  the  numerator  gives  this  slightly 
simpler  form: 


T 


w+pf  Fx(x)dx  +xq  (33) 


£[Z(x)]  = 


pFx(x) 


The  defender  wishes  to  maximize  E[N(T)],  with  x  being  the  decision 
variable.  Therefore  we  must  express  E[N(T)]  in  terms  of  the  expectation  (3.3)  just 
derived. 


12 


If  we  interpret  each  missile  destruction  as  a  point  event,  then  Z(t)  would  be 
the  random  interarrival  time.  If  we  had  a  "continuous  supply"  of  missiles,  and  a 
large  time  window   ,  T,  then  from  the  elementary  renewal  theorem  [Ref.  4]: 


E[N(T)]  - .  (3.4) 

E[Z(x)] 


In  fact,  a  correction  term  might  be  added  [Ref.  4]: 


E[N(T)] 


T       +  Var[Z(z)]-E2[Z(x)]  (3  5) 

E[Z(x))  2£2[Z(t)] 

This  requires  the  computation  of  the  variance  of  Z(x)  which  is  more  difficult  and 
leads  to  a  complex  expression.  In  this  paper  we  retain  the  approximation  (3.4), 
accepting  the  consequent  loss  of  accuracy.  Thus  we  use 


TpFAx) 
E[N(T)]  =  — 


w  +/?  J  Fx(x)dx  +  zq 


(3.6) 


o 


The  defender's  problem  is  to  find: 


:°  =  argmax  E[N(T)] 
0  <  i  s  T    • 


(3.7) 


x°,  the  optimal  threshold  can  be  looked  for  among  the  critical  points  on  the 
interior  of  its  domain,  or  at   x  =  T. 

For  example,  say  X  -  Exponential {X) .   After  substitution  of  the  appropriate 
expressions  for  the  distribution: 
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W)]»  r/(1-gAT>         ;  (3.8) 

X 

Simulations  were  conducted  to  study  the  dependence  of  E[N(T)]  on  the 
parameters  A.,  w,  t,  p.  Plots  resulting  from  these  simulations  are  displayed  and 
discussed  on  Appendix  B. 

Again  we  seek  the  critical  points  of  E[N(T)]  by  finding  the  zeros  of  its 
derivative: 

§ESMA  =0.(i,^*i)«-".i,  (3.9) 

dx  q 

or  equivalently,  by  taking  logarithms  of  both  sides: 

ln(AT+A^  +  i)  =  xx.  (3.10) 

q 

If  p  -  1 ,  it  can  be  verified  that  x°  =  T,  otherwise  the  optimal  threshold  is  to  be 
found  as  a  solution  to  (3.9)  or  (3.10). 

Two  approximate  solutions  for  (3.9)  are  provided.  The  first  is  obtained  by 
rewriting  (3.9)  as 

e^-l  =  At+  —  ,  (3.11) 

q 

and  expanding  the  exponential  term: 

e^  n  i+At+-^.  (3.12) 

2 

Substituting  (3.12)  into  (3.11)  and  solving  for  x   yields: 
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x°  « 


2w 


(3.13) 


M  Xq' 

A  second,  more  accurate  solution,  is  derived  by  taking  (3.13)  as  the  initial 
guess  and  performing  one  iteration  of  Newton's  method  to  solve  (3.1 1),  resulting  in: 


z°  * 


/ 

A 

2wX 

- 1 

/2wA 
V      Q             ^W         1 

e  '    q    + +  1 

(3.14) 


k(e 


2wk 


1) 


Optimal  threshold  values  obtained  with  the  approximations  (3.13)  and  (3.14) 
are  shown  in  Table  3.1,  together  with  numerical  solution  of  (3.10). 


Table  3.1    OPTIMAL  THRESHOLD  APPROXIMATIONS  (p=0.85) 

w\    X 

0.01 

0.05 

0.1 

5.0 

(3.10) 

71.895584 

28.074877 

18.164104 

(3.13) 

81.649658 

36.514837 

25.819889 

(3.14) 

72.721021 

29.928094 

20.659227 

10.0 

(3.10) 

96.921939 

36.328208 

22.991621 

(3.13) 

115.470054 

51.639778 

36.514837 

(3.14) 

99.285331 

41.318454 

29.264013 

20.0 

(3.10) 

128.640050 

45.983242 

28.435654 

(3.13) 

163.299316 

73.029674 

51.639778 

(3.14) 

135.311964 

58.528026 

42.703706 

Although  the  approximations  are  not  always  accurate,  they  are  still  useful  since 
the  expected  number  of  missiles  killed  is  not  very  sensitive  to  the  threshold  value 
(see  Appendix  B).  Table  3.2  shows  values  of    E[N(T)],  obtained  from  (3.8)  with 
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optimal  thresholds  given  by  (3.10),  (3.13)  and  (3.14).  As  it  is  seen,  the  three 
methods  of  deriving  the  optimal  threshold  produced  very  similar  results. 


Table  3.2    EXPECTATION  OF  MISSILES  KILLED  (T=700,  £=0.85) 

w 

\    X 

0.01 

0.05 

0.1 

5.0 

(3.10) 

5.138864 

20.368877 

33.569134 

(3.13) 

5.133375 

20.176869 

32.874955 

(3.14) 

5.138821 

20.357661 

33.476053 

10.0 

(3.10) 

4.777647 

16.784567 

25.374494 

(3.13) 

4.764086 

16.437478 

24.395642 

(3.14) 

4.777394 

16.738027 

25.100826 

20.0 

(3.10) 

4.271479 

12.687247 

17.364481 

(3.13) 

4.241009 

12.197821 

16.343659 

(3.14) 

4.270129 

12.550413 

16.865152 
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C.      SIMULATION 

Optimal  values,  E°[N(T)],  obtained  from  (3.8)  and  (3.10),  were  compared 
to  those  of  simulations,  for  several  values  of  X .  The  simulation  code  is  very 
similar  to  the  program  "Threshold  Simulation"  given  in  Appendix  A.  The  results  are 
displayed  in  Figure  3.1. 


Visible  Kill 


2.5- 


P 
O 
X 
CO 

co 

CO 
CO 

LU 


7.5- 


0.5- 


0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09   0.1 

RATE 


MODEL 


%?%    SIMULATION 


Figure  3.1  Threshold  Policy  -  Results 

As  can  be  seen,  the  model  slightly  underestimates  the  simulated  values, 
although  its  adequacy  is  seen  to  improve  for  larger  values  of  the  killing  rate. 
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IV.   DEFENSE  WITH  PERFECT  TASK  INFORMATION 

A.      PROBLEM  STATEMENT 

Suppose  one  is  given  n  tasks,  with  the  objective  of  completing  the  maximum 
possible  number  of  them  within  a  fixed  time  window  of  length  T.  (E.g.,  a  batch  of 
jobs  arriving  at  a  repair  facility  at  the  beginning  of  the  day;  a  cluster  of  missiles  that 
must  be  destroyed  in  limited  time).  Assume,  as  an  example,  that  the  duration  of  each 
task  is  known,  and  drawn  from  an  Exponential  distribution.  Since  the  objective  is 
to  maximize  the  number  of  tasks  completed,  one  obvious  strategy  is  to  process  the 
tasks  in  ascending  order  of  their  completion  times.  This  policy  will  be  referred  to 
as  a  Smallest-First  (SF)  rule,  in  the  same  fashion  Coffman  [Ref.  5]  uses  Largest- 
First.  In  the  context  of  our  missile  defense  problem  the  above  formulation  is 
artificial  because  task  times  would  never  be  known  in  advance.  Note,  though,  that 
application  of  the  rule  only  requires  knowledge  of  the  order,  and  not  the  actual  task 
(shoot-down)  time.  Therefore  the  results  provide  an  upper  bound  on  the 
effectiveness  of  the  outer  defense  layer,  Dl. 

Denote  by  Xt  the  time  for  completion  of  the  i,h  task.  These  times  are 
independent  and  identically  distributed  exponentially  with  parameter  X  ,  their  i,h 
order  statistic  being  X,~ .  If  the  SF  rule  is  employed,  then  for  some  j  <n: 

*a) +  -  +  xw*  r<  x<i, +  -  +x«  +  w         (4-1} 

in  which  case  only  /'  tasks  may  be  completed  or  "packed"  in  time  T.  Call  the 
number  of  tasks  packed  with  a  SF  rule  N(T),  a  random  variable.  The  aim  of  this 
chapter  is  to  find  the  probability  mass  function  (PMF)  of  N(T),  or  accurate 
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approximations  to  that  PMF.  This  can  then  be  used  to  evaluate  measures  of  system 
effectiveness,  which  can  be  compared  to  the  results  of  applying  other  rules  that  are 
based  on  less  information. 

Similar  problems  have  been  treated  in  recent  literature,  although  in  a  more 
general  way  and  with  different  objectives  and  methodology  than  that  presented  here. 
For  example,  Coffman,  Fayolle,  Jacquet  and  Robert  [Ref.  5]  obtained  asymptotic 
results  for  the  expected  number  packed,  under  a  Largest-First  (LF)  rule;  with  the 
samples  being  drawn  from  a  (7(0,1)  distribution.  Coffman,  Flatto  and  Weber 
[Ref.  6]  derived  the  asymptotic  value  of  the  expected  number  of  intervals  selected 
with  an  optimal  'threshold'  rule,  for  a  class  of  distributions.  Theirs  is  an  "ON- 
LINE" rule,  i.e.,  according  to  their  own  definition,  the  decision  to  pack  or  not  the 
interval  (job)  is  made  at  the  time  of  inspection.  The  missile  problem,  i.e.  of  rank- 
ordering  task  times,  would  also  necessarily  have  to  be  solved  "ON-LINE". 

B.      EXACT  PMF  OF  N(T) 

The  following  well-known  results  are  relevant  to  the  derivation  of  the  PMF 
of  N(T). 

1.       Distribution  of  Spacings 

The  distribution  of  the  differences  between  exponential  order  statistics  is 
next  presented.  More  detail  can  be  found,  e.g.,  in  Barlow  and  Proschan  [Ref.  7]. 

Let  X,,...  yXn  be    i.i.d.    exponential   with    parameter   X.   Defme    the 
spacings,  the  differences  between  successive  order  statistics,  as: 
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Dk  =  X(k)  -  *(*-»  >    with  X(0)  =  °  '  (42) 

k  =  l,...,/i  . 


The  normalized  spacings  nDp  (rt-l)D2, ...  ,Dn  are  again  i.z.d.  exponential  with 
rate  X  .  Let  these  normalized  spacings  be  denoted  by  Y.,...,Y  ,  respectively. 
They  can  be  expressed  in  terms  of  the  order  statistics  in  the  following  way: 

Y{  =  nX{l) 

Y2  =  (n-l)X(2)  -  (n-l)X(1)  (4.3) 

yy  =  (»-;  +  l)Xw  -  (n-;  +  l)XH), 
The  last  equations,  if  solved  for  the  order  statistics,  yield: 

x(d  =  ;ri 

x(2)  --nY^+irj2  (4-4) 

Define  the  summation  of  the  order  statistics: 

'      ft     (°      h  n-i  +  l    *  (4.5) 

j  =  1, ...  ,n. 

Taking  the  expectation  and  variance  of  the  right-hand  side  of  (4.5)  produces: 
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E[S,]  -  ± 


l£,-l-l 


Var[Sj] 


X  ,,i    n-i  +  1 


'j-i  +  1 


(4.6) 


a 


X2i-\  \n-i  +  l) 


2.       Convolution  of  Exponentials 

The  following  result,  for  the  convolution  of  independent  Exponential 
densities,  will  be  used  in  deriving  the  distribution  of  S .   See,  e.g..  Feller  [Ref.  8]. 
Let  Xj , ... ,  X    have  densities 


A,f 


fx(t)  =  X.e  "*  ,  r>0;    kt  +  kki  Vi*k. 
Then  the  sum  Xx  +  ...  +  X    has  density: 


fsW 


(  J        \(  J 


with  the  B.     defined  by  the  following  product: 


*,,,  =  II  <  v K) 

k*i 


(4.7) 


(4.8) 


(4.9) 


3.       Derivation  of  the  PMF  of  N(T) 

In  the  current  case  S  ,  as  defined  in  (4.5),  is  a  linear  combination  of 
j  i.i.d.  exponential  distributions,  as  seen  in  section  III.B.l.  Also  by  a  simple 
transformation  of  scale: 
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Yt  ~  Exponential  I 


n  -i  + 1 


J-i  +  l 


(4.10) 


thus  S  can  be  viewed  as  a  sum  of  j  independent  exponential  random  variables 
with  different  rates.  Then  applying  result  (4.8)  yields  the  distribution,  after  some 
algebra: 


A,w  =  i  zn 


n 


'n-k+1    n-i  +  lV1    ~kt 


k*i 


j  =  1 it-1  ;  A  >  0  ;  t  >  0  . 

To  derive  the  PMF  of  N(T)  the  density  above  is  used: 

P(N(T)=j)  =  Fs   (D  -  F5(D  ; 
witfc  F5(D  =  P(Sj>T)  =  1  -F5(D  , 


(4.11) 


(4.12) 


where  Fs  (t)  is  the  CDF  of  S    and  F5  (f)  is  obtained  from  the  density  (4.11)  by 
integration: 


FsP  =  jfs,W* 


;)E 


A/  n-it+/     n-i  +  n  * 


n 

k*i 


j-k+1     j-i+l 


j-i  +  l       -mkT 


(4.13) 


n-i  +  1 


defined  for  j  =  1,  ...  ,  n-1  . 


Using    (4.12)    we    have    the    desired    PMF    for  /  =  1,..  ,  n-2.  The 
values  /  =0,  n-1,  n   have  not  been  covered  so  far. 
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P(N(T)=n)  =  P(£X.*r)  =  FS(D  ;  (4-14) 


(4.15) 


i-l 


since  Sn   is  a  sum  of  i.i.d.  exponential  Random  Variables,  then: 

Sn  ~  Gamma  ( n,  A )  ,  A  >  0  ; 

-P(N{T)=n)  =  1   -«-*r£^  ; 

If  no  task  is  completed,  the  shortest  time  to  complete  a  task  must  be 
bigger  that  the  time  window,  i.e.: 

P(N(T)=0)  =  FAT)  =  e  nkT  .  (4.16) 

Then  by  a  generalization  of  the  domain  of  definition  in  (4.13), 

0  ifj  *  0 
(4.13)           i//  =  l,...,n-l 

Vf'CW  •,•  '  (4-1?) 

1  ifjzn  +  \ 

With  this  broader  definition  of  F*(T)    (4.12)  holds  for  the   whole  range  of 
/,  (0, 1, ... ,  n),  and  an  exact  PMF  of  N(T)   has  been  derived. 

C.      APPROXIMATION  TO  THE  PMF  OF  N(T) 

To  recapitulate,  an  exact  but  complicated  formula  for  the  PMF  of  N(T)  has 
been  derived  earlier;  see  (4.17).  We  now  desire  an  approximation,  which  might  be 
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easier  to  use.  This  can  be  achieved  via  the  distribution  (CDF)  of  S    in  the  following 

way: 

P(N(T)=j)  =  P(N(T)<j+l)  -  P(N(T)<j) 

(4.18) 
=  Fs  (T)  -  Fs    (T) . 


1.       Gamma  Approximation 

Analysis  of  simulation  data,  e.g.  using  probability  plots,  suggest  that 
the  5  are  approximately  distributed  Gamma  (see  Appendix  C  for  details).  If  the 
approximation  is  good  enough  it  can  be  used  to  compute  values 
of  P(N(T)  =/").  The  values  of  the  Expectation  and  Variance  of  the  5  are 
available  in  (4.6). 

If  5  is  to  be  distributed  Gamma  with  parameters  shape  a.  and 
rate   p  ,  then: 


E[SA  =  -J-    and    Kar[5f]  =  -±  ;  (4.19) 

P,  P/ 


from  which  the  parameter  values  of  the  approximating  distribution  can  be  derived, 
since  the  values  of  the  expectation  and  variance  are  known: 
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a,  = 


Py  = 


'  n-i+l 


Vi- 1 


(4.20) 


xittL,iu-^y 


1=1 


k  -  i  + 1 


1=1 


rc-i  +  1 


It  was  seen  that  this  approximation  produces  good  results  in  computing  the  PMF  of  N(  T ) 
using  (4.18).  Another  approximation  method  is  presented  next. 
2.       Normal  Approximation 

If  we  assume  that  5  can  be  represented  by  a  Normal  distribution  then 
the  resulting  Normal  approximation, with  $  ( )  denoting  the  probability  mass  in  the 
left  tail  of  a  Standard  Normal,  is: 


P(N(T)=j)  =  <J? 


'  T-EiSj] 


* 


f  T-E[Sj+1] 


Jva7[S})      1/^-jJ 


(4.21) 


for  j '  =  1, ... ,  n-\ ;  and  with  E[S}],  Var[S}]    defined  by  (4.6). 

D.      RESULTS 

Once  more  analytical  results  were  checked  through  a  Monte  Carlo  simulation. 
The  simulation  program  replicated  independent  "packing  experiments".  The  values 
of  "Number  Packed"  (or  tasks  completed)  thus  obtained,  were  then  used  in  an  APL 
function  to  compute  the  "experimental"  distribution  of  N{T).  This  distribution  is 
plotted  together  with  the  corresponding  approximations  and  exact  PMF,  (4.12).  The 
inclusion  of  a  continuity  correction  in  the  Normal  Approximation  provided  no 
substantial  improvement.  The  Normal  approximation  was  seen  to  work  better  than 
portrayed  in  Figure  4.1,  for  cases  with  bigger  values  of  n,  as  should  be  expected. 
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The  Gamma  approximation  produced  results  quite  similar  to  those  of  the  exact 
formula. 


PMF  of  N(T) 


0.35 


c       0.3 

o 

to 

o.    0.25 

E 

S        0.2  i 

o 

2    0.15- 


3       0.1 

o 

£     0.05 

0 


n  =  l0;  T=70 


rate  =  0.05 


3       4       5       6       7 

Number  of  Jobs 


NORMAL  |     |  SIMUL      §§3  EXACT      f§§§    GAMMA 


Figure  4.1  Simulation  and  Analytical  Results 

Suppose  an  unsophisticated  policy  is  used,  so  jobs  are  taken  randomly,  without 
ordering.  In  that  case  N(T)  ~  Poisson(X)  ,  and  E[N(T)]  =  XT.  The  PMF  of 
the  number  of  jobs  completed  in  this  way,  was  compared  with  the  corresponding 
PMF  under  a  SF  rule.  This  gives  a  measure  of  the  benefits  of  employing  the  latter 
policy.  Results  are  summarized  in  Table  4.1.  It  can  be  seen  in  the  bottom  row  that, 
for  this  particular  set  of  parameters,  the  SF  rule  accounts  for  an  increase  in  Expected 
value  of  about  85%. 
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Table  4.1    COMPARISON  OF  RESULTS 
FOR  THE  DISTRIBUTION  OF  THE  NUMBER  OF  TASKS  COMPLETED 


j 

SIMULAT 

NORMAL 

GAMMA 

EXACT 

POISSON 

0 

0 

- 

- 

0 

0.030197 

1 

0 

0 

0 

0 

0.105691 

2 

0 

0 

0.000031 

0.000056 

0.184959 

3 

0.002800 

0.000034 

0.002531 

0.002934 

0.215785 

4 

0.037800 

0.019231 

0.035266 

0.035456 

0.188812 

5 

0.152300 

0.195998 

0.158881 

0.156441 

0.132169 

6 

0.315000 

0.342860 

0.310115 

0.309978 

0.077098 

7 

0.301300 

0.252402 

0.301318 

0.303980 

0.038549 

8 

0.152500 

0.122254 

0.151571 

0.151470 

0.016866 

9 

0.035700 

0.047305 

0.036936 

0.036370 

0.006559 

10 

0.002600 

- 

- 

0.003315 

0.002296 

AVG 

6.487500 

6.284765 

6.458033 

6.493038 

3.500000 

E.      APPLICATION  TO  MISSILE  DEFENSE 

We  now  proceed  to  incorporate  the  results  obtained  from  the  packing  problem 
to  a  missile  attack  situation.  As  before,  the  defender  must  acquire  and  shoot  the 
maximum  number  possible  of  n  incoming  missiles,  within  a  restricted  time  window. 

1.       Assumptions 

The  missiles  are  acquired  with  probability  p,  acquisitions  being 
independent  Bernoulli  trials.  A  fixed  time,  w,  is  spent  for  each  acquisition.  The 
defender,  knowing  the  ordering  of  the  "missile  killing  times",  will  select  them  for 
acquisition  with  a  SF  rule.  It  is  not  recognized  if  acquisitions  are  successful  or  not. 
Consequently,  the  first  time  the  defender  fails  to  acquire  a  missile  the  residual  time 
window  is  expended  and  no  more  missiles  are  shot.  Presumably  the  performance 
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could  be  improved  by  introducing  a  threshold,  as  was  done  earlier,  but  this  option 
has  not  been  studied. 
2.       Model 

Let  N(T)  now  represent  the  number  of  missiles  killed  in  time  T,  under 
the  above  assumptions;  then, 

P(N(T)  */)  =  pjFs(T-jw);  (4.22) 

resulting  in  the  PMF: 

P(N(T)  =;)  =  pjFs  (T-jw)  -pJ+iFs    (T-(j+\)w) .  (4.23) 

Recalling  that  F,  (f)  =  1  -Fv  (f)  can  be  obtained  from  (4.13). 

The  result  (4.23)  applied  to  a  given  set  of  parameters  is  plotted  in  Figure 
4.2  together  with  corresponding  values  obtained  through  simulation. 

For  the  purpose  of  comparison  the  results  of  the  model  were  checked 
against  those  obtained  from  employing  the  unsophisticated  strategy,  i.e.,  taking  the 
missiles  in  random  order.  In  this  case  the  distribution  of  N(T)  is  obtained  in  the 
same  fashion,  by  using  (4.23),  but  with  the  difference  that  now  5  is  distributed 
Gamma  with  shape  j  and  rate  X.  The  model  is  seen  to  agree  reasonably  well 
with  the  "experimental  simulation  data"  as  further  detailed  in  Table  4.2.  It  is  also 
verified  a  substantial  improvement  by  using  the  SF  strategy. 

A  noticeable  feature  of  the  distribution  is  that  it  is  no  longer  unimodal,  unlike 
similar  distributions  for  the  packing  problem.  In  plots  generated  with  different  rate 
values  one  local  minimum  appeared  between  two  maxima,  one  of  them  located  at 
zero.  Figure  4.3  depicts  a  situation  in  which  the  rate  is  0.1  and  20  missiles  are 
present,  all  other  parameters  remaining  equal  to  the  previous  situation.  As  can  be 
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seen  the  minimum  is  now  more  sharply  defmed,  and  the  absolute  maximum  is  at 
zero. 


0.25 


Perfect  Task  Information 

PMFof  N(T) 


01       23456789      10 

Number  of  Missiles 


MODEL  ^   SIMULATION 


Figure  4.2  Model  and  Simulation  Results 
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Table  4.2   PERFECT  TASK  INFORMATION  PMF 


j 

■  ■—               —  ■ ...  i  .-—-■■  . 

SIMULATION 

MODEL 

NO  STRATEGY 

0 

0.148574 

0.150000 

0.182953 

1 

0.117059 

0.127500 

0.238861 

2 

0.113057 

0.108782 

0.259735 

3 

0.111056 

0.106964 

0.191905 

4 

0.180590 

0.171499 

0.091931 

5 

0.205103 

0.212791 

0.028368 

6 

0.107554 

0.104298 

0.005541 

7 

0.016508 

0.017283 

0.000659 

8 

0.000500 

0.000872 

0.000045 

9 

0 

0.000011 

0.000002 

10 

0 

0 

0 

AVG 

3.189096 

3.169751 

1.881847 

30 


3  4   5   6    7  8   9  1011121314151617181920 

Number  of  Missiles 


MODEL  g£2   SIMULATION 


Figure    4 . 3 
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V.   CONCLUSIONS 

A.      MODEL  VERIFICATION 

In  the  preceding  chapters  we  developed  three  models  for  defense  against 
missile  attacks.  For  each  model  analytic  results  were  compared  with  simulations,  for 
some  range  or  combination  of  parameters.  To  summarize,  our  three  models 
correspond  to  the  following  situations,  in  ascending  order  of  information  availability: 

Invisible  Kill 

Unknown:    If  acquisition  successful. 

If  and  when  missiles  shot. 
Strategy:      Allocate  time  evenly  between  all  engaged  missiles. 
Decision:     How  many  missiles  to  engage. 

Visible  Kill 

Unknown:  If  acquisition  successful. 

Known:  Time  when  missile  killed  (if  killed). 

Strategy:  Allocate  a  max  threshold  time  to  each  missile  engagement. 

Decision:  Length  of  threshold. 

Perfect  Task  Information 

Unknown:    If  acquisition  successful. 

Known:        Times  of  destruction  of  each  missile. 

Strategy:      Engage  missiles  with  SF  rule. 

Throughout  our  work,  whenever  numerical  results  were  necessary,  assumptions 

on  the  values  of  the  parameters  had  to  be  taken.  We  have  used  what  seemed  to  be 

"reasonable"  values  of  killing  rates,  acquisition  times,  time  window,  etc.  These, 

however,  may  differ  from  actual  parameter  values  to  be  found  in  existing  weapon 

systems.  For  this  reason  we  now  seek  to  "validate"  the  models  in  more  general 

conditions.  We  wish  to  verify  and  compare  our  models  in  a  reasonable  number  of 

parameter  (factors)  combination  without  too  much  simulation  effort.  The  response 
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variable  to  consider  for  validation/comparison  purposes  will  be  the  expected  number 
of  missiles  destroyed.  The  models  depend  on  five  input  parameters  or 
factors:  n,p,  w,  X,  T.  The  same  distribution  for  the  Killing  Time  will  be  used  -  an 
Exponential.  In  order  to  generate  input  parameter  values  we  chose  a  two  level 
fractional  factorial  design,  2V  ,  [Ref.  9],  which  seems  to  attain  the  above  stated 
objectives.  This  design  corresponds  to  the  input  data  set  displayed  on  Table  5.1. 


Table  5.1    INPUT  DATA 


RUN 

X 

D 

P 

w 

T 

1 

0.05 

50 

1.0 

15.0 

70.0 

2 

0.1 

10 

1.0 

5.0 

200.0 

3 

0.1 

50 

1.0 

15.0 

200.0 

4 

0.05 

10 

0.85 

5.0 

200.0 

5 

0.1 

10 

0.85 

5.0 

70.0 

6 

0.1 

50 

0.85 

15.0 

70.0 

7 

0.05 

10 

1.0 

5.0 

70.0 

8 

0.1 

10 

1.0 

15.0 

70.0 

9 

0.05 

10 

0.85 

15.0 

70.0 

10 

0.1 

10 

0.85 

15.0 

200.0 

11 

0.05 

10 

1.0 

15.0 

200.0 

12 

0.05 

50 

0.85 

5.0 

70.0 

13 

0.05 

50 

1.0 

5.0 

200.0 

14 

0.1 

50 

0.85 

5.0 

200.0 

15 

0.1 

50 

1.0 

5.0 

70.0 

16 

0.05 

50 

0.85 

15.0 

200.0 
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1.       Invisible  Kill 

The  simulation  was  carried  out  in  the  following  way:  For  each  run,  the 
optimal  value  of  the  number  of  missiles  to  engage  was  computed  according  to  the 
model  equation  (2.11).  In  each  run  the  simulation  replicates  enough  independent 
attack  experiments  to  get  estimates  of  the  average  number  of  missiles  killed  and 
95%  confidence  bounds.  The  same  average  was  computed  from  the  model  equation 
(2.13).  Results  are  compared  in  Table  5.2. 

Table  5.2    INVISIBLE  KILL  VERIFICATION 


RUN 

LOWER 

AVERAGE 

UPPER 

MODEL 

1 

1.248880 

1.260090 

1.271298 

1.264240 

2 

7.711665 

7.736362 

7.761059 

7.768698 

3 

2.981409 

3.005450 

3.029490 

2.997817 

4 

4.461138 

4.487008 

4.512877 

4.484883 

5 

2.301197 

2.321421 

2.341643 

2.341143 

6 

0.746798 

0.759386 

0.771993 

0.752077 

7 

1.819565 

1.835803 

1.852040 

1.858953 

8 

0.867102 

0.880358 

0.893614 

0.884797 

9 

1.060356 

1.071748 

1.083138 

1.074604 

10 

3.318501 

3.342890 

3.367278 

3.344489 

11 

3.099832 

3.122495 

3.145157 

3.147755 

12 

1.544501 

1.560270 

1.576038 

1.580111 

13 

5.153650 

5.183533 

5.213414 

5.199824 

14 

6.658032 

6.693637 

6.729240 

6.688979 

15 

2.723549 

2.744591 

2.765632 

2.757285 

16 

2.660858 

2.682709 

2.704558 

2.675591 
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As  further  depicted  in  Figure  5.1,  all  runs  produced  satisfying  results. 
Considering  that  (2.13)  is  an  approximation,  the  model  values  seem  to  be  in 
accordance  with  those  of  the  simulation. 


Figure  5.1   Invisible  Kill  Verification 


2.      Visible  Kill  -  Threshold  Strategy 

We  first  determined  the  optimal  threshold  for  each  run  by  numerical 
solution  of  (3.10),  (when  p-\,q-0  the  optimal  threshold  will  be  T).  The 
simulation  is  executed  for  each  combination  of  factors  (run)  and  optimal  threshold. 
Analytical  values  (model)  are  obtained  from  (3.8).  As  can  be  seen  in  the  following 
table  and  figure,  the  results  obtained  from  the  model  are  not  always  accurate  as 
compared  to  the  simulation.  In  fact  the  model  predicts  an  expectation  of  about  13 
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missiles  destroyed,  when  only  10  are  present  in  the  attacks  corresponding  to  run  2. 
This  is  an  expected  result  considering  the  crude  approximations  taken  in  Chapter  III. 
Otherwise  we  find  the  agreement  model/simulation  to  be  quite  acceptable. 

Table  5.3    VISIBLE  KILL  VERIFICATION 


RUN 

LOWER 

AVERAGE 

UPPER 

MODEL 

1 

2.132568 

2.148101 

2.163632 

1.973662 

2 

9.907267 

9.924894 

9.942521 

13.333333 

3 

7.941874 

7.966734 

7.991593 

7.999999 

4 

5.966601 

6.009499 

6.052395 

5.819679 

5 

3.540552 

3.572495 

3.604436 

3.356913 

6 

1.996812 

2.016756 

2.036670 

2.057663 

7 

3.296708 

3.323905 

3.351101 

2.782671 

8 

2.754568 

2.767375 

2.780181 

2.798467 

9 

1.479404 

1.500816 

1.522227 

1.441650 

10 

5.823129 

5.858012 

5.892895 

5.879036 

11 

5.904632 

5.930957 

5.957281 

5.714175 

12 

2.212835 

2.241029 

2.269221 

2.036888 

13 

8.536298 

8.581685 

8.627071 

7.999927 

14 

9.724030 

9.779253 

9.834475 

9.591181 

15 

5.016041 

5.043329 

5.070617 

4.665247 

16 

4.205375 

4.238362 

4.271349 

4.119000 
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Figure  5.2   Visible  Kill  Verification 

3.       Perfect  Task  Information 

Proceeding  in  an  analogous  way  we  compare  the  expected  number  of 
missiles  shot  obtained  from  simulation,  with  model  results  computed  with: 


E[N{T)]-'£jP<kN{T)-j); 

j-0 


(5.1) 


with  P(N(T)  -j)  being  given  in  (4.23). 

The  fit  between  model  and  simulation  is  seen  to  be  excellent  in  this  case. 
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Table  5.4    PERFECT  TASK  INFORMATION  VERIFICATION 


RUN 

LOWER 

AVERAGE 

UPPER 

MODEL 

1 

3.972124 

3.978036 

3.983947 

3.981028 

2 

9.888283 

9.905581 

9.922877 

9.927288 

3 

11.701859 

11.723443 

11.745027 

11.743412 

4 

4.248149 

4.312770 

4.377340 

4.294308 

5 

3.632813 

3.682852 

3.732889 

3.672681 

6 

2.672961 

2.703552 

2.734142 

2.708608 

7 

5.097260 

5.116289 

5.135318 

5.110362 

8 

3.470990 

3.481924 

3.492857 

3.484863 

9 

2.161409 

2.185453 

2.209497 

2.198917 

10 

4.157371 

4.219392 

4.281412 

4.211194 

11 

7.053313 

7.071574 

7.089834 

7.088154 

12 

4.354237 

4.421083 

4.487927 

4.388366 

13 

19.723969 

19.771378 

19.818771 

20.385315 

14 

5.409484 

5.520819 

5.632153 

5.554571 

15 

10.613348 

10.635718 

10.658088 

10.652422 

16 

4.589483 

4.664865 

4.740247 

4.691323 

B.      COMPARISON  OF  STRATEGIES 

Given  that  not  all  the  models  produced  equally  reliable  results,  we  chose  to 
compare  the  different  strategies  using  simulation.  Again,  the  input  data  set  of  Table 
5.1  is  used  to  generate  the  values  depicted  in  Figure  5.4. 

Clearly  the  SF  strategy  in  the  Perfect  Task  Information  situation  shows  the  best 
overall  performance.  This  is  not  a  surprising  result,  since  this  strategy  benefits  from 
the  highest  level  of  information.  Also  the  Threshold  strategy  comes  next  in  overall 
performance,  and  Invisible  Kill  produces  the  worst  results.  More  unexpected  is  the 
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Figure  5.3   Perfect  Task  Information  Verification 


fact  that  in  three  cases  (runs  4,  10,  14)  Threshold  wins  from  Perfect  Task  Information, 
which  is  even  outperformed  by  the  modest  Invisible  Kill  in  run  14.  Although  these 
exceptions  do  not  alter  the  conclusion  that  models  with  access  to  more  information 
perform  better  overall,  it  is  somewhat  surprising  that  this  situation  does  not  occur  for 
every  combination  of  parameters.  A  closer  look  at  the  conspicuous  runs  reveals  that 
in  all  of  them  p  =  0.85 ,  T  -  200.  Thus  the  explanation  for  the  underperformance  of 
Perfect  Task  Information  must  lie  in  the  fact  that  acquisition  is  still  invisible  (see 
IV.E.l).  This  implies  that,  given  sufficient  time,  the  defender  eventually  fails  to 
acquire.  Unlike  with  Threshold  and  Invisible  Kill  there  is  no  chance  of  reacquisition. 
However,  we  see  that  in  more  favorable  conditions,  particularly  when  acquisition 
occurs  with  probability  one,  the  SF  rule  can  enhance  performance  up  to  five  times 
(run  13). 
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Perfect  Info 


Threshold     -^-    Invisible  Kill 


Figure  5.4   Comparison  of  Strategies 

The  shortcomings  in  the  Perfect  Task  Information  strategy,  can  be  overcome 
by  refining  the  model  assumptions  in  IV.E.l.  The  proposed  refinement  consists  on 
including  a  "threshold"  time,  such  that,  after  this  time  has  elapsed  without  target 
destruction  being  observed,  the  defender  reattempts  acquisition  of  the  same  target. 
Various  thresholds  can  be  used  (e.g.,  the  largest  task  time).  The  most  natural 
thresholds  will  be  the  known  task  times  (in  this  case,  knowledge  only  of  the 
ordering  is  not  sufficient). 

With  this  improvement  the  Perfect  Task  Information  model  clearly  dominates 
for  all  simulated  runs,  as  Figure  5.5  shows.  We  do  not  provide  analytical  results  for 
this  improved  model. 


40 


Translated  into  a  context  where  the  weapons  of  the  defender  consist  on 
interceptor  missiles,  these  findings  seem  to  strongly  favor  a  system  that  enables  the 
defender  to  assign  each  interceptor  according  to  the  shortest  time  for  interception 
(SF). 


25 


Improved  Pert.  Info  -*■-  Threshold 


*-    Invisible  Kill 


Figure  5.5   Comparison  with  Improved  Perfect  Information 
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C.      APPLICATION  EXAMPLE 

Suppose  one  is  planning  a  weapon  system  to  face  a  threat  that  allows  a  time 
window  7\  and  that  we  can  forecast  that  n  missiles  will  be  used  in  each  attack. 
We  have  chosen  one  of  the  three  models  to  be  the  most  adequate  to  describe  the 
weapon  system  that  can  be  built.  Then  we  have,  from  our  analytic  work,  the 
function: 

E[N(T)]  =/a,p,w;r,n).  <5-2) 

The  arguments  T  and  n  are  dictated  by  the  threat,  while  \,p,  w  can  be 
incorporated  in  the  design  of  the  defense  system,  within  technical  feasibility  and 
resource  constraints,  and  are  thus  subject  to  decision.  Further,  suppose  we  can 
estimate  the  cost  functions  of  increasing  the  probability  of 
acquisition  ,  c  ,  increasing  the  rate  of  kill  ,  ck,  and  decreasing  the  acquisition 
time  ,  cw.  The  decision  for  the  optimal  combination  of  parameter  to  design  into  the 
system  can  be  obtained  by  solving  the  nonlinear  program: 

Max     E[N(T)] 
s.t.    ck(\)+cp(p)+cw{w)  <>R 

0<k±X  (5.3) 


max 


0<p<  1 

0  <  w  z  w 


max 
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D.      RECOMMENDATIONS 

The  models  developed  in  this  thesis  may  provide  a  basis  for  further  work 
leading  to  realistic  models  of  missile  attack  situations.  For  this  objective  to  be 
attained,  "real  world"  data  and  system  specific  information  must  be  used  to  validate 
the  models.  Further  investigation  may  be  needed  to  verify  the  nature  of  the 
distribution  of  the  killing  times,  or  the  robustness  of  the  models  to  distributional 
assumptions. 

A  possible  refinement  would  be  to  drop  the  assumption  of  invisible  acquisition 
which  is  common  to  all  three  models  but  heavily  penalizes  the  model  of  Perfect 
Task  Information. 

It  might  also  be  useful  to  consider  other  objective  functions.  For  example,  one 
may  argue  that  in  a  scenario  where  some  missiles  are  almost  certain  to  leak,  it 
makes  more  sense  to  attempt  to  minimize  the  maximum  number  of  leaking  missiles. 
Under  some  circumstances  it  may  be  desirable  to  maximize  the  probability  of  killing 
all  attacking  missiles;  under  others  it  may  be  useful  overall  to  expend  some 
resources  to  pass  information  about  leaking  missiles  (course,  estimated  time  of 
arrival)  to  a  secondary  defense  level,  D2.  There  are  numerous  other  options  and 
opportunities. 
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APPENDIX  A.    FORTRAN  CODES 

In  this  appendix  we  list  the  most  important  programs  and  routines.  The  codes 
are  written  in  Fortran77  and  were  run  on  an  Amdahl  5990-5000,  at  the  NPS 
Computer  Center. 

The  simulation  programs  shown  were  used  to  verify  the  numerical  results 
obtained  with  the  probability  models.  The  version  of  the  Perfect  Task  Information 
model  presented  is  the  one  without  the  threshold  improvement.  A  program  to 
generate  exact  (analytical)  results  of  the  Perfect  Task  Information  is  also  included. 
Exact  computations  for  the  other  two  models  are  simple  and  can  be  done  on  a 
programmable  calculator. 

All  simulations  terminate  either  with  the  destruction  of  the  last  missile  or  when 
the  window  elapses. 

A.      INVISIBLE  KILL 

The  main  program,  "Mass  Attack  with  Invisible  Kill",  implements  the 
simulation  to  verify  the  model  developed  in  Chapter  II.  The  program  executes 
independent  replications  of  the  simulation  of  a  missile  attack,  to  compute  the 
expected  number  of  missiles  killed.  This  is  done  for  the  16  sets  of  parameters 
displayed  in  Table  5.1,  and  with    m°    found  as  shown  in  Chapter  II. 

Two  event  routines  are  present:  "Acquisition"  and  "Shooting".  The  simulation 

also  requires   the   SIMUTIL  package    [Ref.   2],   LLRANDOMII   [Ref.    3],   and 

subroutine  STAT. 
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A  descriptive  list  of  the  variables  follows: 

•  "ProbAcq"  is  the  probability  of  acquisition  denoted  by    p    in  Chapter  II. 

•  "TimeWindow"  is  the  time  available  to  destroy  the  targets,     T . 

•  "TimeAcq"  is  the  fixed  acquisition  or  setup  time,     w . 

•  "TimetoKill"  is  the  time  allowed  to  shoot  down  a  single  target. 

•  "Rate"  is  the  killing  rate,  and  also  represents  the  rate  parameter  of  an 
exponential  distribution,     X . 

•  "Average",  "Variance"  and  "Interval Width"  are  statistics  of  the  number  of 
missiles  killed  per  attack,  taken  from  the  sample  of  all  the  simulated  attacks. 
These  statistics  are  computed  by  subroutine  STAT. 

•  "Lower"  and  "Upper"  are  the  95%  confidence  bounds  for  the  expectation  of 
the  number  of  missiles  killed. 

•  "Iseed"  is  the  seed  for  the  random  number  generator.  It  is  initialized  through 
a  "Read"  statement  at  the  beginning  of  the  simulation. 

•  "NrMissiles"  is  the  number  of  missiles  in  an  attack,    n. 

•  "m"  is  the  optimal  number  of  missiles  to  engage.  It  is  computed  through 
approximation  (2.1 1),  as  described  in  Chapter  II.  It  is  an  input  parameter  to  the 
simulation. 

•  "Processed"  counts  missiles  as  they  are  generated. 

•  "NrAcquired"  counts  the  number  of  successful  acquisitions. 

•  "MissilesDown"  counts  the  number  of  missiles  shot.  It  is  used  as  an  input  by 
subroutine  "STAT"  to  compute  the  average  number  of  missiles  killed. 

•  The  integer  variables  "i"  and  "k"  are  indices  used  in  "Do"  loops. 

•  "Acquired"  is  a  logical  variable  and  has  value  .True,  when  the  missile  being 
processed  has  been  successfully  acquired. 
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•  "Destroyed"  is  a  logical  variable.  It  has  value  .True,  when  the  missile  being 
processed  is  destroyed. 

•  "Tsample"    is   the   random   time   to   destroy   a   missile   sampled   from   an 
exponential  distribution  with  rate     X  . 

•  "Now"  is  the  current  simulation  time.  It  is  controlled  by  subroutine  SIMGO 
[Ref.  2]. 

The  Fortran77  code  for  the  simulation  is  as  follows: 


PROGRAM  Mass  Attack  with  Invisible  Kill 

REAL  ProbAcq,  Time  Window,  TimeAcq,  TimetoKill,  Rate,  Average, 
+Variance,  IntervalWidth,  Lower,  Upper 

INTEGER  Iseed,  NrMissiles,  m,  Processed,  NrAcquired,  k, 
+MissilesDown,  i 

COMMON/  SHARE2/  Processed,  NrAcquired,  MissilesDown 
COMMON/  SHARE3/  Iseed,  ProbAcq,  TimeAcq,  TimetoKill,  Rate, 
+NrMissiles 

COMMON/  SHARE4/  Time  Window 
COMMON/  SHARE5/  Average,  Variance,  IntervalWidth,  Lower,  Upper 

Create  input  and  output  files. 
OPEN(  15,  FILE='/EXPINV  DATA') 
OPEN(  20,  FILE='/EXPINV  OUTPUT') 


* 


* 


uiitialize  the  seed  for  the  random  number  generator 
READ(15,  *)  Iseed 
Do  i  =  1,  16 
Read  one  combination  of  input  parameters  (see  Table  5.1),  and  replicate  the 

*  simulation  for  each  of  them. 

Read  (15,*)  Rate,  NrMissiles,  ProbAcq,  TimeAcq,  TimeWindow,  m 
The  time  assigned  to  shoot  at  each  missile  is  computed  from  the  optimal 

*  number  of  missiles  killed  according  to  equation  (2.2). 
If(  m  .ne.  0  )  then 

TimetoKill  =  TimeWindow  /  m  -  TimeAcq 
Else 

TimetoKill  =  0.0 
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End  If 

Initialize  the  statistics. 
Variance  =  0.0 
Average  =  0.0 
Interval  Width  =  1.0 

Perform  independent  replications  of  the  simulation. 
Do  k  =  1,  15000 

CALL  INIT 

CALL  SIMGO 

CALL  STAT(  MissilesDown,  k) 
End  Do 

Write(20,*)  i,  Lower,  Average,  Upper 
End  Do 

CLOSE(15) 
CLOSE(20) 

Stop 

End 

SUBROUTINE  INIT 


*        Initialize  variables  used  in  the  simulation  before  each  replication. 
REAL  Now,  TimeWindow 

INTEGER  Processed,  NrAcquired,  MissilesDown 

* 

COMMON/  SHARE  1/  Now 

COMMON/  SHARE2/  Processed,  NrAcquired,  MissilesDown 

COMMON/  SHARE4/  TimeWindow 

* 

Now  =  0.0 
Processed  =  0 
NrAcquired  =  0 
MissilesDown  =  0 


CALL  CALANDAR  INITIALIZATION 
Schedule  the  first  acquisition  event 
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CALL  SCHEDULE(  Now,  1) 
* 

*  Schedule  end  of  attack  after  time  window  elapsed 
CALL  SCHEDULE(  Now+TimeWindow,  3) 

Return 

End 

SUBROUTINE  ACQUISITION 

REAL  X,  ProbAcq,  Now,  TimeAcq,  TimetoKill,  Rate 

INTEGER  Iseed,  NrAcquired,  Processed,  MissilesDown,  NrMissiles 

LOGICAL  Acquired 

COMMON/  SHARE  1/  Now 

COMMON/  SHARE2/  Processed,  NrAcquired,  MissilesDown 
COMMON/  SHARE3/  Iseed,  ProbAcq,  TimeAcq,  TimetoKill,  Rate, 
+NrMissiles 

SAVE  SHARE2 

* 

*  Sample  from  a  Uniform(0,l),  to  determine  if  acquisition  was  successful 
CALL  RANNUM(  1,  Iseed,  0.0,  1.0,  0.0,  X) 

Acquired  =  .FALSE. 

If(  X  .le.  ProbAcq  )  Acquired  =  .TRUE. 

If(  Acquired  )  then 

NrAcquired  =  NrAcquired  +  1 

*  Schedule  a  shooting  event  for  the  acquired  missile,  after  acquisition  time 
elapsed. 

CALL  SCHEDULE(  Now+TimeAcq,  2) 
Return 
else 

Processed  =  Processed  +  1 

If(  Processed  .It.  NrMissiles  )  then 

*  If  not  acquired  and  not  all  missiles  have  been  processed,  schedule  another 

*  acquisition. 

CALL  SCHEDULE!  Now+TimeAcq+TimetoKill,  1) 
else 

*  End  the  simulation  if  all  missiles  have  been  processed 

CALL  SCHEDULE(  Now+0.00000000001,  3) 
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End  It- 
End  If 

Return 

End 

SUBROUTINE  SHOOTING 

REAL  Now,  TimeAcq,  ProbAcq,  Rate,  Tsample,  TimetoKill 
INTEGER  MissilesDown,  Processed,  Iseed,  NrAcquired,  NrMissiles 
LOGICAL  Destroyed 

COMMON/  SHARE  1/  Now 

COMMON/  SHARE2/  Processed,  NrAcquired,  MissilesDown 
COMMON/  SHARE3/  Iseed,  ProbAcq,  TimeAcq,  TimetoKill,  Rate, 
+NrMissiles 

SAVE  SHARE2 

Samples  the  task  time  from  an  Exponential  distribution 
CALL  RANNUM(  3,  Iseed,  Rate,  0.0,  0.0,  Tsample) 
Destroyed  =  .FALSE. 

If  the  random  task  time  is  less  than  the  time  assigned  to  shoot  the  missile,  it 
will  be  destroyed. 

If(  Tsample  .le.  TimetoKill  )  Destroyed  =  .TRUE. 
If(  Destroyed  )  MissilesDown  =  MissilesDown  +  1 
Processed  =  Processed  +  1 
If(  Processed  .It.  NrMissiles  )  then 

Schedule  another  acquisition  event 

CALL  SCHEDULE(  Now+TimetoKill,  1) 
else 

if  all  missiles  processed,  end  the  attack 

CALL  SCHEDULE(  Now,  3) 
End  If 

Return 

End 
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B.      THRESHOLD  MODEL 

Most  of  what  was  written  for  the  Invisible  Kill  case  still  applies.  All  variables 
retain  their  meanings.  The  variable  "Threshold"  is  introduced  to  denote  the  optimal 
threshold,  x° ,  which  must  be  supplied  to  the  simulation.  In  the  current  case 
optimal  thresholds  were  computed  by  numerical  solution  of  (3.10). 

Fortran  listing  follows: 

PROGRAM  Threshold  Simulation 

REAL  ProbAcq,  Time  Window,  TimeAcq,  Threshold,  Rate 
+,Average,  IntervalWidth,  Variance,  Lower,  Upper 
INTEGER  Iseed,  NrMissiles,  Replications,  MissilesDown,  k 

COMMON/  SHARE2/  MissilesDown 

COMMON/  SHARE3/  Iseed,  ProbAcq,  TimeAcq,  Threshold,  Rate, 
+NrMissiles 

COMMON/  SHARE4/  TimeWindow 
COMMON/  SHARE5/  Average,  Variance,  IntervalWidth,  Upper,  Lower 

OPEN(  15,  FILE=,A,HRESH  DATA') 
OPEN(  20,  FILE='/TKRESH  OUTPUT') 

Read(15,*)  Iseed 
Do  k  =  1,  16 

*  Read  input  parameters  and  optimal  threshold 

Read(15,*)  Rate,  NrMissiles,  ProbAcq,  TimeAcq,  TimeWindow, 
+  Threshold 

Variance  =  0.0 
Average  =  0.0 
Replications  =  0 
IntervalWidth  =  1.0 

*  Replicates  at  least  10,000  times  and  until  95%  confidence  interval  is  less  than 

*  or  equal  to  5%  of  the  average. 

Do  While(  IntervalWidth  .gt.  0.05* Average  .or. 
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+  Replications  .le.  10000  ) 

Replications  =  Replications  +  1 
CALL  INIT 
CALL  SIMGO 

CALL  STAT(  MissilesDown,  Replications) 
End  Do 

Write(20,*)  k,  Lower,  Average,  Upper 
End  Do 

CLOSE(15) 
CLOSE(20) 

Stop 

End 

SUBROUTINE  ACQUISITION 

REAL  X,  Now,  ProbAcq,  TimeAcq,  Threshold,  Rate 
INTEGER  Iseed,  MissilesDown,  NrMissiles 
LOGICAL  Acquired 

COMMON/  SHARE  1/  Now 
COMMON/  SHARE2/  MissilesDown 

COMMON/  SHARE3/  Iseed,  ProbAcq,  TimeAcq,  Threshold,  Rate, 
+NrMissiles 

Sample  from  a  Uniform(0,l)  to  determine  outcome  of  acquisition. 
CALL  RANNUM(  1,  Iseed,  0.0,  1.0,  0,  X) 
Acquired  =  .FALSE. 
If(  X  .le.  ProbAcq  )  Acquired  =  .TRUE. 
If(  Acquired  )  then 
When  missile  is  acquired  schedule  a  "Shooting"  event 

CALL  SCHEDULE(  Now+TimeAcq,  2) 
else 

Upon  failure  to  acquire  schedule  another  acquisition  in  time  =  current  time  + 
acquisition  time  +  threshold. 

CALL  SCHEDULE(  Now+TimeAcq+Threshold,  1) 
End  If 

Return 
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End 

* 

SUBROUTINE  SHOOTING 

REAL  Now,  Tsample,  Threshold,  Rate,  ProbAcq,  TimeAcq 
INTEGER  MissilesDown,  Iseed,  NrMissiles 
LOGICAL  Destroyed 

COMMON/  SHARE  1/  Now 
COMMON/  SHARE2/  MissilesDown 

COMMON/  SHARE3/  Iseed,  ProbAcq,  TimeAcq,  Threshold,  Rate, 
+NrMissiles 

SAVE  SHARE2 

* 

*  Sample  the  exponential  destruction  time 

*  CALL  RANNUM(  3,  Iseed,  Rate,  0.0,  0.0,  Tsample) 
Destroyed  =  .FALSE. 

*  When  the  sampled  time  is  less  than  the  threshold  time,  the  missile  is  killed 
If(  Tsample  .le.  Threshold  )  Destroyed  =  .TRUE. 

If(  Destroyed  )  then 

MissilesDown  =  MissilesDown  +  1 

*  If  all  missile  have  been  shot  ... 

If(MissilesDown  .eq.  NrMissiles)  then 

Return 
else 

*  If  some  missiles  remain,  schedule  another  acquisition  immediately 

CALL  SCHEDULE(  Now+Tsample,  1) 
End  If 
else 

*  If  the  missile  is  not  destroyed,  schedule  an  acquisition  in  time  =  current  time 

*  +  threshold  time 

CALL  SCHEDULE(  Now+Threshold,  1) 

End  If 

* 

Return 

* 

End 
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C.      PERFECT  TASK  INFORMATION 

We  present  programs  for  the  simulation,  "Perfect  Task  Info  Simulation",  and 
for  the  probability  model,  "Exact  PMF  for  Perfect  Task  Info".  This  simulation, 
unlike  the  preceding,  has  no  discrete  event  subroutines.  It  represents  the  case 
without  the  threshold  improvement  mentioned  in  Chapter  V.  The  improved  case, 
however,  is  simulated  with  minor  modifications  to  the  code  presented  here.  The 
same  applies  to  the  plots  displayed  in  Chapter  IV,  for  the  PMF  of  N(T).  To  generate 
those  plots,  we  also  used  a  small  APL  function  which  computes  the  frequency  of 
the  data,  adapted  from  Thomson  [Ref.  10]. 

The  main  program  variables  are: 

"E"  is  a  real  array  for  storage  of  the  task  times. 

"Sum"  accumulates  the  expended  time. 

"Time"  is  the  time  window,     T. 

"Rate"  is  the  rate  parameter,     X,  of  the  exponential  task  times. 

"p"  is  the  probability  of  success  in  acquisition. 

"w"  is  the  fixed  setup  or  acquisition  time. 

"S"  is  used  to  store  samples  from  a  Uniform  distribution. 

"k",  "kk",  "kkk",  are  used  as  indices  for  "Do"  loops. 

"n"  is  the  number  of  missiles. 
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The    commented    simulation    program    is    listed    next.    Subroutine    SORT 
implements  an  insertion  sorting  algorithm,  which  closely  follows  the  code  presented 
by  Etter  [Ref.  11].  Subroutine  STAT  is  also  used  by  the  previous  simulations. 
PROGRAM  Perfect  Task  Info  Simulation 

Simulation  for  the  improved  case 


REAL  E(500),  Sum,  Time,  Rate,  p,  w,  S,  Lower,  Average,  Upper, 
+Variance 
INTEGER  Iseed,  Replications,  k,  kk,  kkk,  n 

COMMON/  MAINSTAT/  Lower,  Average,  Upper,  Variance 

OPEN(  10,  FTLE='/PACKGT  DATA') 
OPEN(  15,  FILE='/PACKGT  OUTPUT') 

Read(10,*)  Iseed,  Replications 

Do  kkk  =  1,  16 

Read(10,*)  Rate,  n,  p,  w,  Time 
Average=  0.0 
Variance  =  0.0 

Do  kk  =  1,  Replications 

Sample  the  n  task  times  from  an  Exponential 
Do  k  =  1,  n 

CALL  RANNUM(  3,  Iseed,  Rate,  0.0,  0.0,  E(k)  ) 
End  Do 

Sorts  the  task  times  in  ascending  order  of  their  lengths 
CALL  SORT(  n,  E) 
k  =  0 
Sum  =  0.0 

Do  While(  k  .le.  n  .and.  Sum  .It.  Time  ) 
k  =  k+  1 

Sample  from  a  Uniform  to  determine  if  acquisition  successful 
CALL  RANNUM(  1,  Iseed,  0.0,  1.0,  0.0,  S) 
If(S  .le.  p)  then 

If  acquired  pack  one  more  job 
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Sum  =  Sum  +  w  +  E(k) 
Else 

If  not  acquired,  expend  time  window 
Sum  =  Time  +1.0 
End  If 
End  Do 
Update  the  statistics 
CALL  STAT(k-l,kk) 
End  Do 

Write(15,*)  kkk,  Lower,  Average,  Upper 
End  Do 

CLOSE(IO) 
CLOSE(15) 

Stop 

End 

SUBROUTINE  SORT(  n,  E) 


*  Sorts  the  elements  of  the  array  E  in  increasing  order.  E(l)  =  smallest; 

*  E(n)  =  largest. 


INTEGER  n,  kk,  k 
REAL  E(n),  Stor 
LOGICAL  Done 

Do  kk  =  1,  n-1 

If(  E(kk)  .gt.  E(kk+1)  )  then 
Done  =  .FALSE. 
k  =  kk 

Do  While(  .not.  Done  ) 
Stor  =  E(k) 
E(k)  =  E(k+l) 
E(k+1)  =  Stor 
If(  k  .eq.  1  .or.  E(k)  .gt.  E(k-l))  then 

Done  =  .TRUE, 
else 
k  =  k-l 
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End  If 
End  Do 
End  If 
End  Do 

Return 

End 

SUBROUTINE  STAT(  y,  n) 

Computes  sequential  statistics 

INTEGER  y,  n 

REAL  Average,  Variance,  Oldaverage,  Upper,  Intervalwidth, 
+01dvariance,  Lower,  X 

COMMON/  MAINSTAT/  Lower,  Average,  Upper,  Variance 

SAVE  MAINSTAT 

Oldaverage  =  Average 
Oldvariance  =  Variance 
X  =  Real(y) 

If(  n  .eq.  1  )  then 

Average  =  X 

Variance  =  0.0 

Intervalwidth  =  1.0E20 
End  If 

If(n  .ge.  2)  Average  =  (REAL(n-l)*01daverage+X)/Real(n) 

If(n  .eq.  2)  Variance  =  Oldaverage**2+(X**2)-(2.0*(Average**2)) 

If(n  .ge.  3)  then 

Variance  =  (REAL(n-2)/REAL(n-l))  *  Oldvariance  + 
+     01daverage**2  -  (REAL(n)/REAL(n-l))  *  Average**2  + 
+     (X**2)/REAL(n-l) 

End  If 

If(n  .GE.  2)  Intervalwidth  =  2.0*1.96*SQRT(Variance)/SQRT(REAL(n)) 

Lower  =  Average  -  Intervalwidth/2.0 

Upper  =  Average  +  Intervalwidth/2.0 
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Return 

End 


The  program  "Exact  PMF  for  Perfect  Task  Info"  implements  the  analytical 
results  derived  in  Chapter  IV.  It  should  be  run  with  option  "(NOXUFLOW",  to 
avoid  underflow  messages. 

Additional  variables  in  this  program  are: 

•  "PMF"  -  represents    P(N(T)=j). 

•  "FSj"  -denotes    PiS^T). 

The  remaining  variables  retain  their  meanings  and  notation. 

PROGRAM  Exact  PMF  for  Perfect  Task  Info 

* 

DOUBLE  PRECISION  Rate,  Time,  p,  w,  PMF,  FSj,  Average 

INTEGER  j,  n,  k 

* 

Open(  20,  File  =  '/PEX  DATA') 
Open(  25,  File  =  '/PEX  OUTPUT') 

Dok=  1,  16 

Read(20,*)  Rate,  n,  p,  w,  Time 
Average  =  0.0D0 
Do  j  =  0,  n 

*  Uses  equation  (4.23)  to  compute  PMF 

PMF  =  (p**j)*(1.0D00-FSj(j,Time-DBLE(j)*w,n,Rate))  - 
+         (p**(j+l))*(1.0D00-FSj(j+l,Time-DBLE(j+l)*w,n,Rate)) 

*  Average  is  computed  with  (5.1) 

Average  =  Average  +  DBLE(j)  *  PMF 
End  Do 

Write(25,*)  k,  Average 
End  Do 
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Close(20) 
Close(25) 

Stop 

End 

DOUBLE  PRECISION  FUNCTION  FSj(jj,X,nn,Rrate) 

Implements  equation  (4.17) 

INTEGER  jj,  i,  k 

DOUBLE  PRECISION  X,  Sum,  Prod,  Rrate,  COMB,  FACT 

If(  X  .le.  0  )  then 
FSj  =  1.0D00 
Return 
End  If 

If(  jj  .le.  0)  then 
FSj  =  0.0D00 
Elself((jj.ge.l)  .and.  (jj.le.(nn-l)))  then 

Sum  =  0.0D00 
Compute  the  product  term  in  equation  (4.13) 
Do  i  =  1,  jj 
Prod  =  1.0D00 
If  ( jj  .ne.  1  )  then 
Do  k  =  1,  jj 

If  (  i  .ne.  k  )  then 

Prod=Prod/(DBLE(nn-k+ 1  )/DBLE(jj-k+ 1  )-DBLE(nn-i+ 1 )/ 
+  DBLE(jj-i+l)) 

End  If 
End  Do 
End  If 

If(ABS(Rrate*(DBLE(nn-i+l)/DBLE(jj-i+l))*X)  .It.  170.0)  then 
Sum=Sum+Prod*(DBLE(jj-i+ 1  )/DBLE(nn-i+ 1  ))*DEXP(-Rrate* 
(DBLE(nn-i+ 1  )/DBLE(jj-i+ 1  ))*X) 
End  If 

FSj  =  Sum  *  COMB(nnjj) 
If(  ABS(FSj)  .gt.  1.0D00)  FSj  =  1.0D00 
End  Do 
Elself(jj  .eq.  nn)  then 
Sum  =  1.0D00 
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Do  k  =  1,  nn-1 

Sum  =  Sum  +  ((Rrate*X)**k)/FACT(k) 

End  Do 

FSj  =  Sum  *  DEXP(-Rrate*X) 
Elself(jj  .gt.  nn)  then 

FSj  =  1.0D00 
End  If 

Return 

End 

DOUBLE  PRECISION  FUNCTION  COMB(  m,  p) 


*  Calculates  the  binomial  coefficients 

*  \P 

INTEGER  m,  p,  k 

COMB  =  1.0D+00 

If(  p  .ne.  m)  then 
Do  k  =  1,  p 

COMB  =  COMB*(DBLE(m-k+l)/DBLE(k)) 
End  Do 

End  If 

* 

Return 

End 

* 

DOUBLE  PRECISION  FUNCTION  FACT(X) 

* 

*  Returns  the  value  of  X! 
INTEGER  X,  K 

FACT  =  1.0D+00 
If(  X  .ge.  2  )  then 

Do  K  =  2,  X 

FACT  =  FACT  *  DBLE(K) 

End  Do 
End  If 
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Return 
End 
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APPENDIX  B.    THRESHOLD  MODEL  GRAPHICS 

This  appendix  contains  plots  for  the  Threshold  model.  The  response  surfaces 
shown  were  obtained    in  GRAFSTAT  [Ref.  12],  from  simulation  output. 

In  both  plots  the  most  salient  feature  is  that  the  expected  number  of  missiles 
killed  is  not  very  sensitive  to  the  threshold  alone. 
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Figure  Bl.   Threshold  and  Acquisition  Time  as  factors 
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THRESHOLD 


p  =  0.75;  n  =  50; 
T  =  500;  w  =  10 
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Figure  B2.   Rate  and  Threshold  as  factors 
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APPENDIX  C.    DISTRIBUTION  OF  Sj 


In  this  appendix,  we  describe  the  evidence  which  suggested  that  a  Gamma 
distribution  would  provide  a  suitable  approximation  to  the  distribution  of  Sj  (the 
approximation  is  developed  in  subsection  IV.C.l). 

A  slight  adaptation  of  the  program  to  simulate  the  packing  problem  (Perfect 
Task  Information)  was  used,  to  produce  samples  of  Sj  for  several  values  of  j.  These 
samples  of  size  1000,  were  used  to  construct  Gamma  probability  plots.  One  such 
plot  is  depicted  in  Figure  CI,  for  j=5. 

Other  plots,  with  different  values  of  j,  produced  much  similar  results.  As  can 
be  seen,  the  sampled  distribution  seems  to  closely  follow  the  shape  of  a  Gamma.  In 
all  plots  slight  departures  from  linearity  were  only  observed  at  the  right  tail. 
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GAMMA  PROBABILITY  PLOT,  N=1000 
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Figure    CI.       Probability    Plot 
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